This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)
## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.
library(tidyr)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(ggplot2)
library(forecast)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(fable)
## Loading required package: fabletools
data <- read.csv("daily_weather.csv")
summary(data)
## number air_pressure_9am air_temp_9am avg_wind_direction_9am
## Min. : 0.0 Min. :908.0 Min. :36.75 Min. : 15.50
## 1st Qu.: 273.5 1st Qu.:916.5 1st Qu.:57.28 1st Qu.: 65.97
## Median : 547.0 Median :918.9 Median :65.72 Median :166.00
## Mean : 547.0 Mean :918.9 Mean :64.93 Mean :142.24
## 3rd Qu.: 820.5 3rd Qu.:921.2 3rd Qu.:73.45 3rd Qu.:191.00
## Max. :1094.0 Max. :929.3 Max. :98.91 Max. :343.40
## NA's :3 NA's :5 NA's :4
## avg_wind_speed_9am max_wind_direction_9am max_wind_speed_9am
## Min. : 0.6935 Min. : 28.90 Min. : 1.186
## 1st Qu.: 2.2488 1st Qu.: 76.55 1st Qu.: 3.067
## Median : 3.8713 Median :177.30 Median : 4.944
## Mean : 5.5083 Mean :148.95 Mean : 7.020
## 3rd Qu.: 7.3372 3rd Qu.:201.23 3rd Qu.: 8.948
## Max. :23.5550 Max. :312.20 Max. :29.841
## NA's :3 NA's :3 NA's :4
## rain_accumulation_9am rain_duration_9am relative_humidity_9am
## Min. : 0.0000 Min. : 0.0 Min. : 6.09
## 1st Qu.: 0.0000 1st Qu.: 0.0 1st Qu.:15.09
## Median : 0.0000 Median : 0.0 Median :23.18
## Mean : 0.2031 Mean : 294.1 Mean :34.24
## 3rd Qu.: 0.0000 3rd Qu.: 0.0 3rd Qu.:45.40
## Max. :24.0200 Max. :17704.0 Max. :92.62
## NA's :6 NA's :3
## relative_humidity_3pm
## Min. : 5.30
## 1st Qu.:17.39
## Median :24.38
## Mean :35.34
## 3rd Qu.:52.06
## Max. :92.25
##
data <- na.omit(data)
air_pressure_9am <- data$air_pressure_9am
air_temp_9am <- data$air_temp_9am
avg_wind_direction_9am <- data$avg_wind_direction_9am
avg_wind_speed_9am <- data$avg_wind_speed_9am
max_wind_direction_9am <- data$max_wind_direction_9am
max_wind_speed_9am <- data$max_wind_speed_9am
rain_accumulation_9am <- data$rain_accumulation_9am
rain_duration_9am <- data$rain_duration_9am
relative_humidity_9am <- data$relative_humidity_9am
relative_humidity_3pm <- data$relative_humidity_3pm
N <- length(relative_humidity_3pm)
#Standardizing the data and converting to time series object
air_pressure_9am <- air_pressure_9am - mean(air_pressure_9am)
air_pressure_9am <- ts(air_pressure_9am, start = 1, end = N)
air_temp_9am <- air_temp_9am - mean(air_temp_9am)
air_temp_9am <- ts(air_temp_9am, start = 1, end = N)
avg_wind_direction_9am <- avg_wind_direction_9am - mean(avg_wind_direction_9am)
avg_wind_direction_9am <- ts(avg_wind_direction_9am, start = 1, end = N)
avg_wind_speed_9am <- avg_wind_speed_9am - mean(avg_wind_speed_9am)
avg_wind_speed_9am <- ts(avg_wind_speed_9am, start = 1, end = N)
max_wind_direction_9am <- max_wind_direction_9am - mean(max_wind_direction_9am)
max_wind_direction_9am <- ts(max_wind_direction_9am, start = 1, end = N)
max_wind_speed_9am <- max_wind_speed_9am - mean(max_wind_speed_9am)
max_wind_speed_9am <- ts(max_wind_speed_9am, start = 1, end = N)
rain_accumulation_9am <- rain_accumulation_9am - mean(rain_accumulation_9am)
rain_accumulation_9am <- ts(rain_accumulation_9am, start = 1, end = N)
rain_duration_9am <- rain_duration_9am - mean(rain_duration_9am)
rain_duration_9am <- ts(rain_duration_9am, start = 1, end = N)
relative_humidity_9am <- relative_humidity_9am - mean(relative_humidity_9am)
relative_humidity_9am <- ts(relative_humidity_9am, start = 1, end = N)
relative_humidity_3pm <- relative_humidity_3pm - mean(relative_humidity_3pm)
relative_humidity_3pm <- ts(relative_humidity_3pm, start = 1, end = N)
plot(air_pressure_9am, xlab="days", ylab = "air_pressure", main = "Air Pressure at 9AM")
plot(air_temp_9am, xlab="days", ylab = "air_temperature", main = "Air Temperature at 9AM")
plot(avg_wind_direction_9am, xlab="days", ylab = "avg_wind_direction_9am", main = "Average Wind Direction at 9AM")
plot(avg_wind_speed_9am, xlab="days", ylab = "avg_wind_speed_9am", main = "Average Wind Speed at 9AM")
plot(max_wind_direction_9am, xlab="days", ylab = "max_wind_direction_9am", main = "Max wind direction at 9AM")
plot(max_wind_speed_9am, xlab="days", ylab = "max_wind_speed_9am", main = "Max Wind Speed at 9AM")
plot(rain_accumulation_9am, xlab="days", ylab = "rain_accumulation_9am", main = "Rain Accumulation at 9AM")
plot(rain_duration_9am, xlab="days", ylab = "rain_duration_9am", main = "Rain Duration at 9AM")
plot(relative_humidity_9am, xlab="days", ylab = "relative_humidity_9am", main = "Relative Humidity at 9AM")
plot(relative_humidity_3pm, xlab="days", ylab = "relative_humidity_3pm", main = "Relative Humidity at 3PM")
acf(air_pressure_9am)
acf(air_temp_9am)
acf(avg_wind_direction_9am)
acf(avg_wind_speed_9am)
acf(max_wind_direction_9am)
acf(max_wind_speed_9am)
acf(rain_accumulation_9am)
acf(rain_duration_9am)
acf(relative_humidity_9am)
acf(relative_humidity_3pm)
#Differencing the variables to make it stationary
diff_air_pressure_9am <- diff(air_pressure_9am)
plot(diff_air_pressure_9am, xlab="days", ylab = "air_pressure_9am", main = "Differenced Air Pressure at 9am Plot")
diff_air_temp_9am <- diff(air_temp_9am)
plot(diff_air_temp_9am, xlab="days", ylab = "air_temp_9am", main = "Differenced Air Temperature at 9am Plot")
diff_avg_wind_direction_9am <- diff(avg_wind_direction_9am)
plot(diff_avg_wind_direction_9am, xlab="days", ylab = "avg_wind_direction_9am", main = "Differenced Average Wind Direction at 9am Plot")
diff_avg_wind_speed_9am <- diff(avg_wind_speed_9am)
plot(diff_avg_wind_speed_9am, xlab="days", ylab = "avg_wind_speed_9am", main = "Differenced Average Wind Speed at 9am Plot")
diff_max_wind_direction_9am <- diff(max_wind_direction_9am)
plot(diff_max_wind_direction_9am, xlab="days", ylab = "max_wind_direction_9am", main = "Differenced Mean Reading at Useless Load Plot")
diff_max_wind_speed_9am <- diff(max_wind_speed_9am)
plot(diff_max_wind_speed_9am, xlab="days", ylab = "max_wind_speed_9am", main = "Differenced Max Wind Speed at 9am Plot")
diff_rain_accumulation_9am <- diff(rain_accumulation_9am)
plot(diff_rain_accumulation_9am, xlab="days", ylab = "rain_accumulation_9am", main = "Differenced Rain Accumulation at 9am Plot")
diff_rain_duration_9am <- diff(rain_duration_9am)
plot(diff_rain_duration_9am, xlab="days", ylab = "rain_duration_9am", main = "Differenced Rain Duration at 9am Plot")
diff_relative_humidity_9am <- diff(relative_humidity_9am)
plot(diff_relative_humidity_9am, xlab="days", ylab = "relative_humidity_9am", main = "Differenced Relative Humidity at 9am Plot")
diff_relative_humidity_3pm <- diff(relative_humidity_3pm)
plot(diff_relative_humidity_3pm, xlab="days", ylab = "relative_humidity_3pm", main = "Differenced Relative Humidity at 3pm Plot")
adf.test(diff_air_pressure_9am)
## Warning in adf.test(diff_air_pressure_9am): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_air_pressure_9am
## Dickey-Fuller = -16.148, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(diff_air_temp_9am)
## Warning in adf.test(diff_air_temp_9am): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_air_temp_9am
## Dickey-Fuller = -15.332, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(diff_avg_wind_direction_9am)
## Warning in adf.test(diff_avg_wind_direction_9am): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_avg_wind_direction_9am
## Dickey-Fuller = -17.135, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(diff_avg_wind_speed_9am)
## Warning in adf.test(diff_avg_wind_speed_9am): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_avg_wind_speed_9am
## Dickey-Fuller = -15.752, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(diff_max_wind_direction_9am)
## Warning in adf.test(diff_max_wind_direction_9am): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_max_wind_direction_9am
## Dickey-Fuller = -16.261, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(diff_max_wind_speed_9am)
## Warning in adf.test(diff_max_wind_speed_9am): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_max_wind_speed_9am
## Dickey-Fuller = -15.698, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(diff_rain_accumulation_9am)
## Warning in adf.test(diff_rain_accumulation_9am): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_rain_accumulation_9am
## Dickey-Fuller = -16.832, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(diff_rain_duration_9am)
## Warning in adf.test(diff_rain_duration_9am): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_rain_duration_9am
## Dickey-Fuller = -16.862, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(diff_relative_humidity_9am)
## Warning in adf.test(diff_relative_humidity_9am): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_relative_humidity_9am
## Dickey-Fuller = -16.118, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(diff_relative_humidity_3pm)
## Warning in adf.test(diff_relative_humidity_3pm): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_relative_humidity_3pm
## Dickey-Fuller = -15.81, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
#ACF
acf(diff_air_pressure_9am)
acf(diff_air_temp_9am)
acf(diff_avg_wind_direction_9am)
acf(diff_avg_wind_speed_9am)
acf(diff_max_wind_direction_9am)
acf(diff_max_wind_speed_9am)
acf(diff_rain_accumulation_9am)
acf(diff_rain_duration_9am)
acf(diff_relative_humidity_9am)
acf(diff_relative_humidity_3pm)
#PACF
pacf(diff_air_pressure_9am)
pacf(diff_air_temp_9am)
pacf(diff_avg_wind_direction_9am)
pacf(diff_avg_wind_speed_9am)
pacf(diff_max_wind_direction_9am)
pacf(diff_max_wind_speed_9am)
pacf(diff_rain_accumulation_9am)
pacf(diff_rain_duration_9am)
pacf(diff_relative_humidity_9am)
pacf(diff_relative_humidity_3pm)
#CCF
ccf(diff_relative_humidity_3pm, diff_air_pressure_9am)
ccf(diff_relative_humidity_3pm, diff_air_temp_9am)
ccf(diff_relative_humidity_3pm, diff_avg_wind_direction_9am)
ccf(diff_relative_humidity_3pm, diff_avg_wind_speed_9am)
ccf(diff_relative_humidity_3pm, diff_max_wind_direction_9am)
ccf(diff_relative_humidity_3pm, diff_max_wind_speed_9am)
ccf(diff_relative_humidity_3pm, diff_rain_accumulation_9am)
ccf(diff_relative_humidity_3pm, diff_rain_duration_9am)
ccf(diff_relative_humidity_3pm, diff_relative_humidity_9am)
spectrum(diff_air_pressure_9am, log='no', xlab='Frequencies',
ylab='Power',
main='Periodogram Air Pressure at 9am')
spectrum(diff_air_temp_9am, log='no', xlab='Frequencies',
ylab='Power',
main='Periodogram Air Temperature at 9am')
spectrum(diff_avg_wind_direction_9am, log='no', xlab='Frequencies',
ylab='Power',
main='Periodogram Average Wind Direction at 9am')
spectrum(diff_avg_wind_speed_9am, log='no', xlab='Frequencies',
ylab='Power',
main='Periodogram Average Wind Speed at 9am')
spectrum(diff_max_wind_direction_9am, log='no', xlab='Frequencies',
ylab='Power',
main='Periodogram Max Wind Direction at 9am')
spectrum(diff_max_wind_speed_9am, log='no', xlab='Frequencies',
ylab='Power',
main='Periodogram Max Wind Speed at 9am')
spectrum(diff_rain_accumulation_9am, log='no', xlab='Frequencies',
ylab='Power',
main='Periodogram Rain Accumulation at 9am')
spectrum(diff_rain_duration_9am, log='no', xlab='Frequencies',
ylab='Power',
main='Periodogram Rain Duration at 9am')
spectrum(diff_relative_humidity_9am, log='no', xlab='Frequencies',
ylab='Power',
main='Periodogram Relative Humidity at 9am')
spectrum(diff_relative_humidity_3pm, log='no', xlab='Frequencies',
ylab='Power',
main='Periodogram Relative Humidity at 3pm')
suppressWarnings({
#Based on observations from ACF and CCF creating new lagged inputs and output
air_pressure_9am_lag1 <- filter(diff_air_pressure_9am, c(0,1), sides=1)
air_pressure_9am_lag2 <- filter(air_pressure_9am_lag1, c(0,1), sides=1)
air_temp_9am_lag1 <- filter(diff_air_temp_9am, c(0,1), sides=1)
air_temp_9am_lag2 <- filter(air_temp_9am_lag1, c(0,1), sides=1)
avg_wind_direction_9am_lag1 <- filter(diff_avg_wind_direction_9am, c(0,1), sides=1)
avg_wind_direction_9am_lag2 <- filter(avg_wind_direction_9am_lag1, c(0,1), sides=1)
avg_wind_speed_9am_lag1 <- filter(diff_avg_wind_speed_9am, c(0,1), sides=1)
avg_wind_speed_9am_lag2 <- filter(avg_wind_speed_9am_lag1, c(0,1), sides=1)
max_wind_direction_9am_lag1 <- filter(diff_max_wind_direction_9am, c(0,1), sides=1)
max_wind_direction_9am_lag2 <- filter(max_wind_direction_9am_lag1, c(0,1), sides=1)
max_wind_speed_9am_lag1 <- filter(diff_max_wind_speed_9am, c(0,1), sides=1)
max_wind_speed_9am_lag2 <- filter(max_wind_speed_9am_lag1, c(0,1), sides=1)
rain_accumulation_9am_lag1 <- filter(diff_rain_accumulation_9am, c(0,1), sides=1)
rain_accumulation_9am_lag2 <- filter(rain_accumulation_9am_lag1, c(0,1), sides=1)
rain_duration_9am_lag1 <- filter(diff_rain_duration_9am, c(0,1), sides=1)
rain_duration_9am_lag2 <- filter(rain_duration_9am_lag1, c(0,1), sides=1)
relative_humidity_9am_lag1 <- filter(diff_relative_humidity_9am, c(0,1), sides=1)
relative_humidity_9am_lag2 <- filter(relative_humidity_9am_lag1, c(0,1), sides=1)
relative_humidity_3pm_lag1 <- filter(diff_relative_humidity_3pm, c(0,1), sides=1)
relative_humidity_3pm_lag2 <- filter(relative_humidity_3pm_lag1, c(0,1), sides=1)
diff_air_pressure_9am <- diff_air_pressure_9am[1:N]
air_pressure_9am_lag1 <- air_pressure_9am_lag1[1:N]
air_pressure_9am_lag2 <- air_pressure_9am_lag2[1:N]
diff_air_temp_9am <- diff_air_temp_9am[1:N]
air_temp_9am_lag1 <- air_temp_9am_lag1[1:N]
air_temp_9am_lag2 <- air_temp_9am_lag2[1:N]
diff_avg_wind_direction_9am <- diff_avg_wind_direction_9am[1:N]
avg_wind_direction_9am_lag1 <- avg_wind_direction_9am_lag1[1:N]
avg_wind_direction_9am_lag2 <- avg_wind_direction_9am_lag2[1:N]
diff_avg_wind_speed_9am <- diff_avg_wind_speed_9am[1:N]
avg_wind_speed_9am_lag1 <- avg_wind_speed_9am_lag1[1:N]
avg_wind_speed_9am_lag2 <- avg_wind_speed_9am_lag2[1:N]
diff_max_wind_direction_9am <- diff_max_wind_direction_9am[1:N]
max_wind_direction_9am_lag1 <- max_wind_direction_9am_lag1[1:N]
max_wind_direction_9am_lag2 <- max_wind_direction_9am_lag2[1:N]
diff_max_wind_speed_9am <- diff_max_wind_speed_9am[1:N]
max_wind_speed_9am_lag1 <- max_wind_speed_9am_lag1[1:N]
max_wind_speed_9am_lag2 <- max_wind_speed_9am_lag2[1:N]
diff_rain_accumulation_9am <- diff_rain_accumulation_9am[1:N]
rain_accumulation_9am_lag1 <- rain_accumulation_9am_lag1[1:N]
rain_accumulation_9am_lag2 <- rain_accumulation_9am_lag2[1:N]
diff_rain_duration_9am <- diff_rain_duration_9am[1:N]
rain_duration_9am_lag1 <- rain_duration_9am_lag1[1:N]
rain_duration_9am_lag2 <- rain_duration_9am_lag2[1:N]
diff_relative_humidity_9am <- diff_relative_humidity_9am[1:N]
relative_humidity_9am_lag1 <- relative_humidity_9am_lag1[1:N]
relative_humidity_9am_lag2 <- relative_humidity_9am_lag2[1:N]
diff_relative_humidity_3pm <- diff_relative_humidity_3pm[1:N]
relative_humidity_3pm_lag1 <- relative_humidity_3pm_lag1[1:N]
relative_humidity_3pm_lag2 <- relative_humidity_3pm_lag2[1:N]
time_prd <- 1:N-1
sin_pred <- sin(2 * pi * time_prd * (1/365))
cos_pred <- cos(2 * pi * time_prd * (1/365))
data_arima2 <- cbind(air_pressure_9am_lag1, air_pressure_9am_lag2,
air_temp_9am_lag1, air_temp_9am_lag2,
avg_wind_direction_9am_lag1, avg_wind_direction_9am_lag2,
avg_wind_speed_9am_lag1, avg_wind_speed_9am_lag2,
rain_duration_9am_lag1, rain_duration_9am_lag2,
relative_humidity_9am_lag1, relative_humidity_9am_lag2,
relative_humidity_3pm_lag1, relative_humidity_3pm_lag2, sin_pred, cos_pred)
data_arima1 <- cbind(air_pressure_9am_lag1, air_temp_9am_lag1, avg_wind_direction_9am_lag1, avg_wind_speed_9am_lag1, rain_duration_9am_lag1, relative_humidity_9am_lag1, relative_humidity_3pm_lag1, sin_pred, cos_pred)
})
for(i in 1:5) {
for (j in 1:5){
complex_model_lag2 <- arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,], order=c(i,0,j))
pr <- predict(complex_model_lag2, newxreg = data_arima2[1001:1050,])
ts.plot(diff_relative_humidity_3pm[1000:1049], pr$pred, lty=1:2, col=c("blue", "red"))
legend("topleft", legend = c("Actual", "Predicted"), col = 1:2, lty = 1:2)
cat()
cat("AIC:", AIC(complex_model_lag2), "BIC:", BIC(complex_model_lag2), "p:", i, "q:", j,"\n")
cat()
}
}
## AIC: 9065.332 BIC: 9163.447 p: 1 q: 1
## AIC: 9065.844 BIC: 9168.865 p: 1 q: 2
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,
## : possible convergence problem: optim gave code = 1
## AIC: 9064.37 BIC: 9172.296 p: 1 q: 3
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,
## : possible convergence problem: optim gave code = 1
## AIC: 9067.453 BIC: 9180.285 p: 1 q: 4
## AIC: 9070.327 BIC: 9188.065 p: 1 q: 5
## AIC: 9066.041 BIC: 9169.062 p: 2 q: 1
## AIC: 9067.595 BIC: 9175.521 p: 2 q: 2
## AIC: 9066.132 BIC: 9178.964 p: 2 q: 3
## AIC: 9065.692 BIC: 9183.431 p: 2 q: 4
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,
## : possible convergence problem: optim gave code = 1
## AIC: 9066.021 BIC: 9188.665 p: 2 q: 5
## AIC: 9067.933 BIC: 9175.86 p: 3 q: 1
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## AIC: 9068.679 BIC: 9181.512 p: 3 q: 2
## AIC: 9065.597 BIC: 9183.335 p: 3 q: 3
## AIC: 9067.818 BIC: 9190.462 p: 3 q: 4
## AIC: 9068.575 BIC: 9196.124 p: 3 q: 5
## AIC: 9069.878 BIC: 9182.711 p: 4 q: 1
## Warning in log(s2): NaNs produced
## AIC: 9070.332 BIC: 9188.07 p: 4 q: 2
## AIC: 9068.148 BIC: 9190.792 p: 4 q: 3
## AIC: 9069.889 BIC: 9197.439 p: 4 q: 4
## AIC: 9060.747 BIC: 9193.203 p: 4 q: 5
## AIC: 9071.647 BIC: 9189.385 p: 5 q: 1
## AIC: 9071.674 BIC: 9194.318 p: 5 q: 2
## AIC: 9069.755 BIC: 9197.304 p: 5 q: 3
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,
## : possible convergence problem: optim gave code = 1
## AIC: 9064.824 BIC: 9197.279 p: 5 q: 4
## AIC: 9069.728 BIC: 9207.09 p: 5 q: 5
for(i in 1:5) {
for (j in 1:5){
# Fit ARMAX model
armax_model <- arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,], order=c(i, 0, j), include.mean = FALSE)
# Predict using ARMAX model
pr_armax <- predict(armax_model, newxreg = data_arima1[1001:1050,])
# Plot actual vs predicted values
ts.plot(diff_relative_humidity_3pm[1000:1049], pr_armax$pred, lty=1:2, col=c("blue", "red"))
legend("topleft", legend = c("Actual", "Predicted"), col = 1:2, lty = 1:2)
cat()
# Print AIC and BIC scores
cat("AIC:", AIC(armax_model), "BIC:", BIC(armax_model), "p:", i, "q:", j,"\n")
cat()
}
}
## AIC: 9069.574 BIC: 9128.455 p: 1 q: 1
## AIC: 9071.266 BIC: 9135.054 p: 1 q: 2
## AIC: 9073.089 BIC: 9141.783 p: 1 q: 3
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## AIC: 9074.745 BIC: 9148.346 p: 1 q: 4
## AIC: 9075.388 BIC: 9153.896 p: 1 q: 5
## AIC: 9071.029 BIC: 9134.817 p: 2 q: 1
## AIC: 9073.199 BIC: 9141.894 p: 2 q: 2
## AIC: 9074.595 BIC: 9148.196 p: 2 q: 3
## AIC: 9070.652 BIC: 9149.16 p: 2 q: 4
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1
## AIC: 9072.456 BIC: 9155.87 p: 2 q: 5
## AIC: 9072.848 BIC: 9141.542 p: 3 q: 1
## AIC: 9073.779 BIC: 9147.38 p: 3 q: 2
## AIC: 9075.753 BIC: 9154.261 p: 3 q: 3
## AIC: 9072.472 BIC: 9155.887 p: 3 q: 4
## AIC: 9072.756 BIC: 9161.078 p: 3 q: 5
## AIC: 9074.385 BIC: 9147.986 p: 4 q: 1
## AIC: 9075.912 BIC: 9154.42 p: 4 q: 2
## AIC: 9077.974 BIC: 9161.389 p: 4 q: 3
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1
## AIC: 9073.77 BIC: 9162.092 p: 4 q: 4
## AIC: 9073.691 BIC: 9166.919 p: 4 q: 5
## AIC: 9075.776 BIC: 9154.284 p: 5 q: 1
## AIC: 9077.228 BIC: 9160.643 p: 5 q: 2
## AIC: 9074.709 BIC: 9163.03 p: 5 q: 3
## AIC: 9076.222 BIC: 9169.45 p: 5 q: 4
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1
## AIC: 9069.652 BIC: 9167.787 p: 5 q: 5
for(i in 1:5) {
for (j in 1:5){
complex_model_lag1 <- arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,], order=c(i,0,j))
pr <- predict(complex_model_lag1, newxreg = data_arima1[1001:1050,])
ts.plot(diff_relative_humidity_3pm[1000:1049], pr$pred, lty=1:2, col=c("blue", "red"))
legend("topleft", legend = c("Actual", "Predicted"), col = 1:2, lty = 1:2)
cat()
cat("AIC:", AIC(complex_model_lag1), "BIC:", BIC(complex_model_lag1), "p:", i, "q:", j,"\n")
cat()
}
}
## AIC: 9070.219 BIC: 9134.007 p: 1 q: 1
## AIC: 9071.847 BIC: 9140.541 p: 1 q: 2
## AIC: 9073.646 BIC: 9147.247 p: 1 q: 3
## AIC: 9075.057 BIC: 9153.566 p: 1 q: 4
## AIC: 9075.86 BIC: 9159.275 p: 1 q: 5
## AIC: 9071.572 BIC: 9140.267 p: 2 q: 1
## AIC: 9073.509 BIC: 9147.11 p: 2 q: 2
## AIC: 9075.181 BIC: 9153.689 p: 2 q: 3
## AIC: 9071.24 BIC: 9154.655 p: 2 q: 4
## AIC: 9073.071 BIC: 9161.393 p: 2 q: 5
## AIC: 9073.323 BIC: 9146.924 p: 3 q: 1
## AIC: 9074.454 BIC: 9152.962 p: 3 q: 2
## AIC: 9076.22 BIC: 9159.635 p: 3 q: 3
## AIC: 9072.999 BIC: 9161.32 p: 3 q: 4
## AIC: 9071.952 BIC: 9165.18 p: 3 q: 5
## AIC: 9075.012 BIC: 9153.52 p: 4 q: 1
## AIC: 9076.385 BIC: 9159.8 p: 4 q: 2
## AIC: 9078.449 BIC: 9166.771 p: 4 q: 3
## AIC: 9074.308 BIC: 9167.536 p: 4 q: 4
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1
## AIC: 9071.216 BIC: 9169.351 p: 4 q: 5
## AIC: 9076.413 BIC: 9159.828 p: 5 q: 1
## AIC: 9077.747 BIC: 9166.069 p: 5 q: 2
## AIC: 9075.244 BIC: 9168.472 p: 5 q: 3
## AIC: 9076.712 BIC: 9174.847 p: 5 q: 4
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1
## AIC: 9069.659 BIC: 9172.701 p: 5 q: 5
for(i in 1:5) {
for (j in 1:5){
simple_model <- arima(diff_relative_humidity_3pm[0:1000], order=c(i,0,j))
cat()
cat("AIC:", AIC(simple_model), "BIC:", BIC(simple_model), "p:", i, "q:", j,"\n")
cat()
}
}
## AIC: 9074.406 BIC: 9094.037 p: 1 q: 1
## AIC: 9076.126 BIC: 9100.665 p: 1 q: 2
## AIC: 9077.978 BIC: 9107.425 p: 1 q: 3
## AIC: 9079.405 BIC: 9113.759 p: 1 q: 4
## AIC: 9079.756 BIC: 9119.018 p: 1 q: 5
## AIC: 9075.475 BIC: 9100.014 p: 2 q: 1
## AIC: 9077.663 BIC: 9107.109 p: 2 q: 2
## AIC: 9073.938 BIC: 9108.293 p: 2 q: 3
## AIC: 9068.856 BIC: 9108.118 p: 2 q: 4
## AIC: 9069.232 BIC: 9113.402 p: 2 q: 5
## AIC: 9076.505 BIC: 9105.952 p: 3 q: 1
## AIC: 9078.539 BIC: 9112.893 p: 3 q: 2
## AIC: 9074.797 BIC: 9114.059 p: 3 q: 3
## AIC: 9080.866 BIC: 9125.036 p: 3 q: 4
## Warning in arima(diff_relative_humidity_3pm[0:1000], order = c(i, 0, j)):
## possible convergence problem: optim gave code = 1
## AIC: 9080.726 BIC: 9129.803 p: 3 q: 5
## Warning in log(s2): NaNs produced
## AIC: 9078.304 BIC: 9112.658 p: 4 q: 1
## AIC: 9079.957 BIC: 9119.219 p: 4 q: 2
## AIC: 9076.99 BIC: 9121.159 p: 4 q: 3
## AIC: 9078.994 BIC: 9128.071 p: 4 q: 4
## AIC: 9078.601 BIC: 9132.587 p: 4 q: 5
## AIC: 9079.83 BIC: 9119.092 p: 5 q: 1
## AIC: 9081.592 BIC: 9125.761 p: 5 q: 2
## AIC: 9078.971 BIC: 9128.049 p: 5 q: 3
## AIC: 9080.99 BIC: 9134.975 p: 5 q: 4
## AIC: 9081.844 BIC: 9140.737 p: 5 q: 5
complex_model_lag2 <- arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,], order=c(5, 0, 4))
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima2[0:1000,
## : possible convergence problem: optim gave code = 1
average_model_lag1 <- arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,], order=c(5, 0, 5))
## Warning in arima(diff_relative_humidity_3pm[0:1000], xreg = data_arima1[0:1000,
## : possible convergence problem: optim gave code = 1
simple_model <- arima(diff_relative_humidity_3pm[0:1000], order=c(2, 0, 4))
AIC(complex_model_lag2)
## [1] 9064.824
BIC(complex_model_lag2)
## [1] 9197.279
AIC(average_model_lag1)
## [1] 9069.659
BIC(average_model_lag1)
## [1] 9172.701
AIC(simple_model)
## [1] 9068.856
BIC(simple_model)
## [1] 9108.118
#Residuals are independent
complex_resids <- residuals(complex_model_lag2)
Box.test(complex_resids,lag=3)
##
## Box-Pierce test
##
## data: complex_resids
## X-squared = 0.0013766, df = 3, p-value = 1
plot(complex_resids, ylab='Residuals',main='Residual plot for complex model 2')
#Residuals are independent
complex_resids <- residuals(complex_model_lag2)
Box.test(complex_resids,lag=3)
##
## Box-Pierce test
##
## data: complex_resids
## X-squared = 0.0013766, df = 3, p-value = 1
plot(complex_resids, ylab='Residuals',main='Residual plot for complex model 2')
pr <- predict(complex_model_lag2, newxreg = data_arima2[1001:1050,])
ts.plot(diff_relative_humidity_3pm[1000:1049], pr$pred, lty=1:2, col=c("blue", "red"),
ylab='relative_humidity_3pm', main='Actual Vs Predicted')
legend("topleft", legend = c("Actual", "Predicted"), col = 1:2, lty = 1:2)
pr <- predict(complex_model_lag2, newxreg = data_arima2[1001:1050,])
ts.plot(diff_relative_humidity_3pm[1000:1049], pr$pred, lty=1:2, col=c("blue", "red"),
ylab='relative_humidity_3pm', main='Actual Vs Predicted')
legend("topleft", legend = c("Actual", "Predicted"), col = 1:2, lty = 1:2)